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Abstract 



We have developed a new numerical code to study the evolution of distorted, 
rotating black holes. We discuss the numerical methods and gauge conditions 
we developed to evolve such spacetimes. The code has been put through a 
series of tests, and we report on (a) results of comparisons with codes designed 
to evolve non-rotating holes, (b) evolution of Kerr spacetimes for which ana- 
lytic properties are known, and (c) the evolution of distorted rotating holes. 
The code accurately reproduces results of the previous NCSA non-rotating 
code and passes convergence tests. New features of the evolution of rotating 
black holes not seen in non-rotating holes are identified. With this code we 
can evolve rotating black holes up to about t = lOOM, depending on the res- 
olution and angular momentum. We also describe a new family of black hole 
initial data sets which represent rotating holes with a wide range of distortion 
parameters, and distorted non-rotating black holes with odd-parity radiation. 
Finally, we study the limiting slices for a maximally sliced rotating black hole 
and find good agreement with theoretical predictions. 

PACS numbers: 04.30.+X, 95.30.Sf, 04.25.Dm 
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I. INTRODUCTION 



In the last few years there has been much study of black holes. They are expected to be 
an important source of gravitational waves. The LIGO and VIRGO interferometers should 
begin taking data later this decade and sensitive bar detectors are already on line and 
will improve 0. At the same time computer power is increasing dramatically, permitting 
well resolved numerical simulations of axisymmetric black holes and the first 3D simulations 
of black holes. Collectively this work will lead ultimately to numerical simulations of 3D 
spiraling, coalescing black holes, which will be essential to interpreting the gravitational 
waveforms expected to be detected. 

Astrophysical black holes in the universe are expected to possess angular momentum, 
yet due to the difficulty of black hole simulations most numerical calculations of vacuum 
black hole initial data sets, such as those discovered by Bowen and York ||^ nearly 15 years 
ago, have not been attempted until now. (Recently, however, rotating matter configurations 
have been successfully collapsed to the point where a horizon forms [Q.) Rotation adds 
a new degree of freedom in the system, complicating matters significantly. Even in the 
stationary case, the Kerr solution is much more complex than the Schwarzschild solution. 
As all black hole systems with net angular momentum must eventually settle down to a 
perturbed Kerr spacetime, it is essential to develop techniques to study distorted rotating 
black holes numerically. 

Such studies are interesting not only because they allow one to examine the nonlinear 
evolution of distorted, rotating single black holes, but also because they should be useful 
in understanding the intermediate and late coalescence phase of the general collision of two 
rotating black holes, as both cases correspond to highly distorted Kerr spacetimes. This 
parallel between distorted single black holes and the collision of two black holes was striking 
in the non- rotating case 0-0], and we expect the same to be true for the rotating case. 

In this series of papers we show how one can construct and evolve vacuum, distorted, 
rotating black holes. In this first paper in this series we present details of a numerical code 
designed to evolve rotating black hole initial data sets, such as Kerr, Bowen and York [Q, 
and a new class of data sets we have developed 0. With such a code we can now study 
the dynamics of highly distorted, rotating black holes, paralleling recent work of the NCSA 
group on non-rotating holes P,pi]|. The first studies of the physics of these evolving systems 



are presented in a companion paper ||ri[, referred to henceforth as Paper II, and a complete 
study of the initial data sets will be discussed in Ref. @], which we refer to as Paper III. 

The paper is organized as follows: In section |I| we review the basic formulation of the 
equations, discuss our choice of variables, and analyze symmetries present in our spacetimes. 
The initial data sets we have used are discussed in section |T|, and the numerical code and 
gauge choices are described in section In section we describe the behavior of the 
metric functions in a rotating code and discuss the differences with the non-rotating code. 



II. PRELIMINARIES 
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A. 3+1 Formalism 



We use the well known ADM [[T^ formulation of the Einstein equations as the basis for 



our numerical code. Pertinent details are summarized here, but we refer the reader to [T^ for 
a complete treatment of this formalism. In the ADM formalism spacetime is foliated into a 
set of non-intersecting spacelike surfaces. There are two kinematic variables which describe 
the evolution between these surfaces: the lapse a, which describes the rate of advance of 
time along a timelike unit vector normal to a surface, and the spacelike shift vector f]'^ 
that describes the motion of coordinates within a surface. The choice of lapse and shift is 
essentially arbitrary, and our choices will be described in section [IV A] below. 
The line element is written as 



- PaP'')dr + 2Padx''dt + -fabdx^dx" 



(1) 



where 'jab is the 3-metric induced on each spacelike slice. Given a choice of lapse a and 
shift vector (3^, the Einstein equations in 3+1 formalism split into evolution equations for 
the 3-metric jab and constraint equations that must be satisfied on any time slice. The 
evolution equations, in vacuum, are 

dtjij = -2aK,j + V,/5, + V, A (2) 
dtKi, = -V.V,« + a (Ri, + K Ki, - 2K,„irf ) + r^mK^j + K,^W,r + K^j^i^, (3) 

where Kij is the extrinsic curvature of the 3-dimensional time slice. The Hamiltonian 
constraint equation is 



R + K^ - K'^Kij = 
and the three momentum constraint equations are 

Vi (k'^ - y^k) = 0. 



(4) 



(5) 



Note that in the above, K = KijY\ Pi = P^'jij, and Vj is the covariant three-space deriva- 



tive. As we discuss in section |T|, the constraint equations are used to obtain the initial 
data, and the evolution equations are used to advance the solution in time. 



B. Definition of Variables 



We build on earlier work of Ref. in defining the variables used in our code. Addi- 



tional metric and extrinsic curvature variables must be introduced to allow for the odd-parity 
modes present now that the system allows for rotation. We define the variables used in our 
evolutions as follows: 



Iriri l-qd ^ri4> 



A C Esm^e 



li, = I 7r,e lee 760 C B FsinO \ (6a) 

l'n<p le<t> l^<t> J \E sin^ 9 FsinO D sin^ 9 



and 
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/ Ha He He sin^ 9 \ 
Kij = <iI^Hij = \ He Hb Hp sine . (6b) 
V He sin^ 9 He sin 9 Hd sin^ 9 J 

In these expressions 1] is a. logarithmic radial coordinate, and {9,(j)) are the usual angular coor- 
dinates. The relation between rj and the standard radial coordinates used for Schwarzschild 
and Kerr black holes is discussed in section |T|. As in Ref. 0, the conformal factor is 
determined on the initial slice and, since we do not use it as a dynamical variable, it re- 
mains fixed in time afterwards. The introduction of \E' into the extrinsic curvature variables 
simplifies the evolution equations somewhat. For the purposes of our numerical evolution 
we will treat A, B, D, F, all six if's, and all three components of the shift as dynamical 
variables to be evolved. Two of the shift components are used to eliminate metric variable 
C (this will be discussed in section IV A 2|) , and one shift component is used to eliminate E. 



The various factors of sin6' are included in the definitions to explicitly account for some of 
the behavior of the metric variables near the axis of symmetry and the equator. 

Within this coordinate system, the ADM mass and angular momentum about the z-axis 



are defined to be 14 



MaDM = „ r 



^ jVa{'^e-'''^)dS^ (7a) 



In terms of the variables defined in this paper these expressions yield 

Madm = - r 6"/^ {dr,^ - ^/2) sin 9d9, (8a) 



J = P^ = ^ r ^^HeJ ^ sin^ 9d9 (8b) 
4 Jo V yi 

Because of this, the variable He is extremely important. It determines whether angular 
momentum is present in the spacetime. Although the ADM mass is defined strictly only 
at spatial infinity I^, in practice we evaluate it at the edge of the spatial grid. As we 
use a logarithmic radial coordinate rj, this is in the asymptotic regime. While the angular 
momentum is, in principle, also measured at J°, the presence of the azimuthal Killing vector 
makes it possible to evaluate J at any radius. We compute this quantity during our evolution 
and use it as a test of the accuracy of our code. This will be discussed in a later section. 



C. Symmetries 

Symmetries are an important consideration in the evolution code for the setting of bound- 
ary conditions. This is important not only for solving numerical elliptic equations, but also 
for appropriate finite differencing of our variables. Most of the appropriate conditions can 
be derived merely by considering the behavior of metric functions near the boundaries or 
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the symmetry operations one can perform on a spinning object, without appeahng to the 
Einstein equations themselves. 
The principal symmetries are: 

(i) Axisymmetry. We chose to study axisymmetric spacetimes as a first step towards 
understanding general rotating black holes. Even in the stationary, Kerr case, rotating black 
holes are already an inherently 2D problem (i.e., they cannot be treated as a spherical system, 
as Schwarzschild is.) The fact that any rotating black hole should settle down to a Kerr black 
hole, and thus an axisymmetric solution, makes this a good choice for understanding the 
late time behavior of any rotating black hole system. One physical restriction imposed by 
this symmetry is that angular momentum cannot be radiated (see, e.g., p. 297 of Ref. [0). 

Axisymmetry requires that the transformation 0^0 + constant leaves the problem 
unchanged, and therefore all variables in the problem must be independent of 0. However, 
it is important to realize that — >■ —0 is not a symmetry of this spacetime, since performing 
this transformation would amount to reversing the direction of spin of the hole. In polar 
coordinates, there is an additional consequence of axisymmetry: the transformation — >■ 
+ 7r produces the same result as 6 ^ —6, providing a boundary condition on the symmetry 
axis ^ = 0. These considerations require the following variables to be symmetric across the 
axis: A, B, D, F Ha, Hb, Hjy, Hp, Hp, a, P^, and P^. The remaining variables C and He 
are antisymmetric across the axis. 

(m) Equatorial Plane Symmetry. We require the spacetime to be identical when reflected 
through the equator defined hj 6 = 7i /2 {z = 0). This symmetry was not strictly necessary, 
but adopting it reduces the complexity of the problem slightly. One might think that the 
transformation z —z (or equivalently 6 ^ n — 6) would result in a hole spinning in the 
opposite direction, but it does not. Inverting through the equator in this fashion is actually 
equivalent to rotating the hole by vr radians about the x axis, then by sending — 0. 
Both of these latter operations clearly result in reversing the sense of the hole's spin, and 
so performing both of them leaves the direction of the hole's spin unchanged. One can 
easily check that the Kerr solution itself is manifestly unchanged by the operation z —>■ —z. 
These considerations require the following variables to be symmetric across the equator: A, 
B, D, Ha, Hb, Hd, He, a, P^, and P'^. The remaining variables C, He, Hp, and F are 
antisymmetric across the equator. 

There is an alternate equatorial boundary condition to consider. In the "cosmic screw" 
(the collision of two black-holes on axis with equal and opposite angular momenta) the sense 
of the rotation of the two holes is changed upon inversion through the equatorial plane. The 
appropriate isometry is thus {9, 0) ^ (vr — 9, —0). Hp and P'^ are now antisymmetric, and 
Hp and F are symmetric. These conditions will be considered in a future paper. 

(iii) Time/rotation symmetry. This is a symmetry of the initial slice. It simplifies the 
initial data for the extrinsic curvature by requiring that all but the values Hp and Hp be 
zero. This symmetry says that the transformation (0, t) (— 0, — t) leaves the problem 
unchanged. 

{iv) Inversion through the throat. Building on previous work 0,|ry], we construct our 
spacetimes to be inversion symmetric through the black hole "throat" . This Einstein- Rosen 
bridge [|T6| construction has two geometrically identical sheets connected smoothly at the 
throat, located at r] = 0. This symmetry requires the metric to be invariant under the 
transformation rj —i]. For the case of the Kerr spacetime, it can be expressed as 
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f ^ (m^ - a^) / (4f) , (9) 

where f is a generalization of the Schwarzschild isotropic radius, a is the Kerr rotation 
parameter, and m is the mass of the Kerr black hole. The variable f is defined as 



^^lEZe. (lOa) 



and is related to the usual Boyer-Lindquist [17| radial coordinate via 



f 1 + — — 1 + — — . (10b) 



2f / V 2f 



Note that in the Kerr spacetime the horizon, located at r = m + \/iv? — a^, is at f = 
yjm? — a? /2 in the f coordinates, or at 77 = 0, just as in previous studies of the Schwarzschild 
spacetime . 

This symmetry is perhaps the most important because it is what makes our spacetime 
a black hole spacetime. When combined with symmetry {iii) above and the differential 
equation for a trapped surface it tells us that 17 = is a trapped surface on the initial 
slice. To impose this symmetry, one requires that metric variables with a single rj index 
will be antisymmetric across the throat, and all others will be symmetric across the throat. 
The extrinsic curvature variables will have the opposite symmetry on the throat as the 
corresponding metric variables when the lapse is antisymmetric and the same symmetry 
when the lapse is symmetric. Thus, the following are symmetric: A, B, D, F, He, He, f3^, 
and P'^. The following are antisymmetric: C, E, Ha, Hb, H^, Hp, (this is the symmetry 
of the if's when a is antisymmetric). 

It is possible to slice the Kerr spacetime with a symmetric lapse across the throat, but 
a different symmetry at the throat needs to be employed for certain extrinsic curvature 
variables. Kerr initial data specifies that He and Hp are symmetric and antisymmetric 
across the throat, respectively; Bowen and York ||^ initial data sets, discussed below, specify 
that He is symmetric across the throat. These conditions are incompatible with an 77 — > —r) 
boundary condition if a symmetric lapse is employed {jij must have the same symmetry as 
aHij for the evolution equations to be consistent). This problem could be removed if we 
use (?7, 0) — i> (— ?7, —0) instead of 77 — > —rj, resulting in consistent evolution equations on 
both sides of the throat. We point out that this condition is consistent with the Kerr initial 
data, and is the generalization of the technique of evolving Schwarzschild with a symmetric 
lapse [TO . 



III. INITIAL DATA 

We can use the constraints to construct initial data. An especially convenient formulation 
of the initial data problem was given by Bowen and York 0. They define a conformal 
spacetime in which 

-fij = (11a) 
K,, = ^~^H,,. (lib) 
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In terms of the conformal extrinsic curvature, we write the Kij as 



Ha He He sin^ 9 \ 
He Hb Hp sine 
y He sin^ 9 Hp sin 9 H^ sin^ 9 j 



(12) 
(13) 



If the initial data set is maximally sliced (tri^' = 0), or if the value of tri^' is held constant, 
then the conformal factor drops out of the momentum constraint equations. This decouples 
the constraints and allows one to solve for the momentum constraints first and then to solve 
the Hamiltonian constraint for the conformal factor. For more detail, see Ref. 

The initial data sets we have developed to study distorted rotating black holes are based 
on the both the "Brill Wave plus Black Hole" solutions of the NCSA group and the 
rotating black holes of Bowen and York |Q. A complete discussion and analysis of these 
rotating data sets will be published in Paper III. Here we provide only the basic construction 
of the initial data sets. Generalizing the conformally flat approach of Bowen and York , 
and following earlier work of Ref. [O], we write the metric with a free function q{r], 9): 



dV 



.2(g-<?o) 



(drf + d9'' 



+ sin 



(14) 



The function go is chosen so that if the function g(?7, 9) is chosen to vanish, we are left 
with the Kerr 3-metric in these coordinates, as shown below. (We note that Kerr is not 
conformally flat.) For other choices of q{f],9) we obtain another spacetime. As we will see, 
appropriate choices of this function (along with appropriate solutions to the momentum 
constraint for the extrinsic curvature terms) can lead to Schwarzschild, the NCSA distorted 
non- rotating black hole (as in Ref. ^^), the Bowen and York rotating black hole [Q, or a 
distorted Bowen and York black hole. 

As in Ref. |19[, the function q, representing the Brill wave, can be chosen somewhat 
arbitrarily, subject to symmetry conditions on the throat, axis and equator, and falloff 
conditions at large radii [§,0- Often the function q will be chosen to have an inversion 
symmetric gaussian part given by: 



q = sm''9qG, 
qc = Qo {e~'+ + e 



(15a) 
(15b) 
(15c) 



This form of the Brill wave will be characterized by several parameters: Qo (its amplitude), 
a (its width), tjq (its coordinate location), and n, specifying its angular dependence, which 
must be positive and even. We note that Eqs. (15) simply provide a convenient way to 
parameterize the initial data sets, and to allow us to easily adjust the "Brill wave" part of 
the initial data. Many other devices are possible. 

In this generalization, we may now interpret the parameter gg as the Brill wave q required 
to make the spacetime conformally flat. A distorted Bowen and York spacetime can be made 
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by setting q = qc + Qo and a distorted Kerr spacetime can be made by setting q = sin" Oqc, 
when appropriate solutions are taken for the momentum constraints. 

As a consequence of this generahzation of the metric, we must now solve both the Hamil- 
tonian and momentum constraints. The Hamiltonian constraint equation is 

5^ ^ \ ^-7.. 2 . 2n rr2\ 



(16) 



and the component of the momentum constraint (the only non-trivial momentum con- 
straint equation for the initial data) is 

drjHE sin^ e + de [Hf sin^ O) = 0. (17) 

Based on these ideas, we consider several solutions to the initial value problem for evo- 
lution in this paper. Complete details of these data sets and extensions of them will be 
discussed in Paper III. The first solution we consider is the Kerr solution, where 

^o = ^?i?/sm^^ (18a) 
= = aie' (18b) 

Hp = -2^0" CL^^r COS e sin^ ^p~^ (18c) 
He = ^o" a ^ sin^ ^ (p^ (^^ - + 2r^ (r^ + a^) ) A"^/ ( 18d) 

= r^ + cos^e, (18e) 

A = r^ -2mr + a^. (18f) 

Here is the value of for the Kerr spacetime, and glf^ denotes a Kerr metric element 



in Boyer-Lindquist coordinates, as given in Ref. IT^]. Although this may not look familiar, 
one may check that with the transformations 

r = r+ cosh^ (^/2) — r_ sinh^ (^7/2) , (19a) 

and 

r± = m ± ^w? — (19b) 



the Kerr solution in Boyer-Lindquist coordinates actually results, as detailed in Paper III. 
Note that this is the same transformation given in Eqs.(|lOD. Note also that if a = then 
go = and we recover the Schwarzschild 3-metric. This metric is now in the form used in 
the previous NCSA black hole studies [|,|,|lO|,|l|,|0|. 



Another solution we study is related to the original solution of Bowen and York 
Following them, we choose 



He 
Hp 



3J 




(20a) 
(20b) 
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where J = am is the total angular momentum of the spacetime. But then we may distort 
it by choosing 

q = qGsm^e + qo (21) 

and solving the Hamiltonian constraint. 

A third solution is the odd-parity distorted non-rotating black hole, which we create by 
choosing 

Hj^ = q'^ ((n' + 1) - (2 + n') sin^ o) sin"'-^' 6 (22a) 
Hp = -drjq'a cos 6 sin"'"^ 6, (22b) 



which has "odd-parity" radiation but no rotation, as one can verify by checking Eq. (|8E|). 
Because He vanishes in the large t] limit, the ADM integral for the angular momentum must 
also vanish (See Eq. (pb|)). Also note that n' must be odd and have a value of at least 3. 

Thus, in some sense this third solution is a distortion from a Schwarzschild black hole, 
but one that is fundamentally different from the "Brill wave plus black hole" data sets of 
Ref . . This data set describes a non-rotating spacetime, parts of which are "twisting and 
untwisting," but without any angular momentum. Some of these data sets will be evolved 
and studied in detail in Paper II. 

For all these data sets we need to solve the nonlinear elliptic Hamiltonian constraint 
equation. We accomplish this by replacing \1' with + 6"$ and expanding in 6"$ to create a 
linear elliptic equation and iteratively solving for 6"$. The linear elliptic equation is solved 
by employing finite differencing and the multigrid elliptic solver described in section |1V A . 
For a more detailed discussion, see Paper III. 



IV. THE EVOLUTION CODE 

In this section we discuss our evolution code, including various gauge choices and com- 
putational details. 



A. Gauge Choices 

1. Lapse 

Maximal slicing is used, with a lapse that is antisymmetric across the throat. Max- 
imal slicing is the condition that the trace of the extrinsic curvature is zero throughout 
the evolution. We note that the Kerr solution in standard form is maximally sliced with 
antisymmetric lapse. Setting tiK = in the evolution for tiK gives 

= - VVaa + aR. (23) 

We find that the code is more stable if we use the energy constraint equation to eliminate 
the Ricci scalar R from Eq. (^) since it contains second derivatives of the metric functions. 
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The behavior of second derivatives can be troublesome near the peaks in metric functions 
that can develop near black holes. The evolution equation for tiK becomes, therefore, 

= -V^Vatt + aK'^Kij. (24) 

This equation is solved numerically on each time slice during the evolution using a multigrid 
solver to be discussed in section IIV B. 



Maximal slicing is an important slicing condition for its singularity avoiding properties. 
As the black hole evolves, constant rj observers are really falling into the hole. Observers 
at the same value of r] but different values of (j) along the equator will, for example, find 
that they are closer together as the evolution proceeds. If we were not careful, our code 
would follow these observers as they fall in until the distance between them becomes zero (as 
they reach the singularity) causing our code to crash. For Schwarzschild and distorted non- 
rotating spacetimes maximal slicing prevents the evolution from reaching the singularity. 
Instead, when the slice approaches the radius r = |M the lapse goes to zero and proper 
time will not advance at that location . 

For rotating black holes the singularity structure inside the horizon is significantly dif- 
ferent. For these reasons one might worry that the singularity avoidance properties could 
be different. The interior region of a rotating black hole has a timelike singularity, Cauchy 



horizons, and an infinite patchwork of "other universes" |2^. However, this problem was 
investigated by Eardley and Smarr |2^ , who showed that maximal slicing does in fact avoid 
the singularity inside the Kerr black hole (in fact, it avoids the Cauchy horizon), although 
in a different manner from Schwarzschild. They estimated the location of a limit surface, 
which was found with more precision by Duncan [E^ who obtained 



3M 



/4, (25) 



where is given in terms of the standard Boyer-Lindquist coordinates. However, we will 
measure "the circumferential radius," not the Boyer-Lindquist r coordinate, so the above 
equation must be converted using the following formula: 



+ a2 + 2a'^m/r. (26) 



This is simply the result of evaluating Vc = \I'^-\AD at ^ = n/2 for the analytic Kerr metric. 

Thus the limit surface depends on the rotation parameter. Note that if the rotation pa- 
rameter a vanishes, we recover the well known limiting surface Vc = 3M/2 for Schwarzschild. 
Actual time slices will approach but do not fall inside this radius interior to the horizon. 
However, we usually use an antisymmetric slicing condition at the throat, causing the lapse 
to vanish there, preventing from evolving on throat. Further out from the throat, however, 
the coordinate system will evolve toward this limiting surface. 

We have tested these ideas by numerically evolving rotating black holes with maximal 
slicing using our code, described in detail below. First we evolve a non-rotating hole, but 
with our new black hole construction described above. In Fig. ^ we plot the limit surface 
using the circumferential radius for an odd-parity distorted non-rotating hole, and we see 
that it locks onto the familiar = |M surface. The parameters describing the initial data 
for this calculation are (Qq, n', r/g, a') = (2.0, 3, 1.0, 1.0). This run is described as run ol in 
Paper II. 
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Next, Fig. shows the computed value of Tc for a rapidly rotating {a/m = .70) black 
hole. The calculation of from the data on the slice t = 60M is given with a solid line, 
and the limit surface predicted by Duncan with a dashed line. In terms of the parameters in 
section |ITl| the initial data for the a/m = .70 spacetime is described by {Qq, tiq, a, J, n) = 
(1.0, 1.0, 1.0, 10.0, 2). This same calculation is labeled run in Paper II, where its evolution 
is analyzed in detail. 

As one can see, maximal slicing does indeed serve to limit the advance of proper time 
inside a rotating black hole, in spite of the fact that the singularity structure is so different 
from that of a non-rotating black hole. The slices do converge to the correct radius, as 
predicted by Duncan, near the horizon (r/ ^ 2.7). Outside this region the radius is, of 
course, much larger. This spacetime is not expected to have exactly the same limiting 
structure as the symmetrically sliced Kerr spacetime, since it has been strongly distorted 
by gravitational waves and since an antisymmetric lapse has been employed (freezing the 
evolution at the throat); hence the variations about Duncan's prediction. We have tested 
this prediction for limit surfaces for a number of rotating black holes and find similar results 
in all cases. We note that this observation could be used to extract the rotation parameter, 
a/m, from the maximally sliced spacetime. We also note that these figures are very useful 
for estimating the location of the apparent horizon, even without an apparent horizon finder, 
since the slices will wrap up close to the limit surface inside, but move quickly away from it 
as one moves out, crossing the horizon. 



2. Shift 

In previous work on non-rotating black holes a shift vector was chosen to make the off- 
diagonal component of the 3-metric 7^0 vanish |T0| (in that system no other off-diagonal 
components were present). This condition was found to be crucial to suppress a numerical 
instability occurring near the axis of symmetry. An elliptic equation for the two components 
of the shift {(3^ and f3^) present in that system provided a smooth shift and stable evolution. 
Following this philosophy we have developed a shift condition that generalizes this approach 
to the rotating case. 

The gauge condition used in our evolutions is C = and E = 0. Our choice has the 
property of reducing to the NCSA gauge used in the non-rotating work as the rotation 
goes to zero, and reducing to the Kerr shift as Qo 0. Since the Kerr shift allows the 
stationary rotating black hole metric to be manifestly time independent, one expects that 
for the dynamical case a similar shift will be helpful. The evolution equations for C and 
E can be used to construct a differential equation for the shifts. Other off-diagonal terms 
may be eliminated through appropriate shift choices, and as with the shift used in the non- 
rotating system, there is additional gauge degree of freedom that has not been exploited. 
We note that the quasi-isotropic shift has been used successfully in a recent study of rotating 
matter collapse 0]. 

Let us now consider how to implement the condition C = and E = 0. First, the 
relevant metric evolution equations are: 

dtC = 0= -2aHc + Ad^(3^ + 59^/5" + F sin 68^(3''' (27a) 
dtE = 0= -2aHE + Ddr^(3'*' + F/^^ sin 6. (27b) 
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These equations can be combined to produce a single equation involving and 

2a I^Hc - ^^H^ = Ad,P' + " ^) ^eP"^- (28) 

We can solve this equation by introducing an auxiliary function Q through the defini- 
tions: 

(3"^ = deQ, (29a) 
(3^ = dr,n, (29b) 



(following Ref. [0), producing an elliptic equation for the function Q: 

2a (hc - ^^He^ = Ad'^n + - ^) (30) 

This equation is then solved by finite differencing using a numerical elliptic equation solver 
discussed in the next section. The solution Q is then differentiated by centered deriva- 
tives to recover the shift components and according to Eqs. (p9D. In practice, these 
shifts remain fairly small during the evolution. Their main function is to suppress the axis 



instability, as noted in Ref. ||T9[ where a similar shift was used. 



Once Q is known, P"^ can be calculated by integrating Eq. (|27b| ): 

P'^= r dr](2aHE - FdriP^/ sine) /D. (31) 

Only one boundary condition needs to be set (the outer boundary condition is most conve- 
nient), and it is generally set equal to the Kerr value. The inner boundary condition, that 
P'^ must be symmetric across the throat, is guaranteed by Eq. ( |27bD . This shift component 
is needed to keep the coordinates from becoming "tangled up" as they are dragged around 
by the rotating hole. Without such a shift the coordinates would rotate, leading to metric 



shear [^. This shift component, P''', is typically larger than P''^ or P . 

This method for obtaining the shift has proved effective, although there are some nu- 
merical difficulties that should be mentioned. First, since P'^ is computed by integrating 
an ODE, errors tend to accumulate as the integration progresses inward. This can cause 
trouble when integrating across the sharp peaks that develop in the metric functions near 
the horizon of a black hole (See section ^ below for discussion of these peaks). On the 
other hand, this occurs in the region where the lapse collapses significantly (near the inner 
portion of the grid), and is not noticeable before the axis instability sets in (see section |^ 
for a discussion of the instability), so it has not been a serious problem. For Kerr initial 
data, the numerically computed shift had a maximum deviation from the analytic solution 
of about 0.1% on a 200 x 55 grid. We show a plot of the shift in Fig. ^ for a Kerr black 
hole with a rotation parameter oi a/m = .676. For dynamic black holes, the shift takes on 
a similar form. 

In addition, each 6 =constant line of integration is independent of the others. As a result 
of this decoupling, there can be fiuctuations in the shift across different angular zones near 
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the axis. This becomes more apparent at late times when the axis instabihty, common in 
most axisymmetric codes, sets in. As noted above this problem is not noticeable until after 
the instability sets in, and so it is not a cause of difficulty. An example of this problem is 
shown in the next section. 



B. Computational Issues and Numerical Issues 

The code was developed using MathTensor, a package that runs under Mathematica, to 
convert the equations to Fortran readable form. Scripts have been written that automatically 
produce Fortran code given symbolic input in MathTensor, so different variable choices and 
gauge conditions can be tested fairly easily. 

The actual evolution code was written in Fortran 77 to run on the Cray Y-MP and 
the Cray C-90. Currently, for a 200 x 55 grid it obtains about 160 MFlops on a single 
processor C-90. About 70% of the time to evolve a spacetime with this code is spent solving 
elliptic equations (two elliptic equations must be solved on each time slice). The multigrid 
linear system solver used by our code for the elliptic equations was provided by Steven 
Schaffer of New Mexico Technical Institute p6|| (UMGS2). We find that although it does 
not achieve the highest code performance in the traditional sense (measured by MFlops), 
it produces a solution in very few iterations, so the time to solution is often less than 
we achieve with "higher performance" iterative solvers [^. UMGS2 is a semi-coarsening 
multigrid code. This differs from full multigrid (described in detail in the above reference) 
by only performing the coarsening along the angular dimension of the grid. This is quite 
useful because the spacetimes vary much less in the angular direction than in the radial 
direction. 

The evolution method that we have chosen is leapfrog. Leapfrog is an explicit evolution 
scheme that requires us to keep the metric variables on two time steps, and the extrinsic cur- 
vature data on two other time steps sandwiched between the metric variables. Schematically 
the evolution looks like this: 

3 1 , , 

lt+M/2 = -It - -It-At (32a) 

It+At = lt + ^t'j {Kt+At/2,lt+At/2) (32b) 

~ 3 1 

Kt+At = - ■^Kt-At/2 (32c) 

Kt+3At/2 = Kt+M/2 + ^tk (kt+At, It+At) (32d) 



Spatial derivatives needed in the above equations were calculated using centered, second 
order finite differencing. This scheme is essentially the same as that discussed in Ref. ||10|| . 

As in most asymmetric black hole codes to date, our code has an axis instability (See 
Ref. [|I^ for a detailed discussion of the axis instability). This instability grows worse as 
rotation increases. An effective strategy to slow the growth of this instabihty is to reduce the 
number of angular zones, thus keeping them farther away from the axis. As mentioned above, 
we have also used the gauge freedom in the equations to eliminate off-diagonal elements in 
the 3-metric that tend to exacerbate this instability. However, at late times, (t ~ 70— lOOM), 
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when peaks in metric functions become large, strong instabilities can develop, causing the 
code to crash. These next two plots show the function B for the run labeled r4, a distorted 
Bowen and York (rotating) black hole. In Fig. ^(a) we see the function at t = lOOM for 
a grid resolution of 150 x 24. In Fig. ^h) we see the same function for t = 70M with a 
grid resolution of 300 x 48 (all angular zones are shown in the figure). The instability is 
clearly visible at earlier times in the higher resolution calculation, while the lower resolution 
calculation is still stable at later times. Note that it is the region in the interior of the black 
hole that develops difficulty (the horizon is located near rj = 3, where the dip in the metric 
function B occurs). 

Because the radial resolution is most crucial, we often perform calculations with reso- 
lutions of 300 X 30 zones, providing both good accuracy and late time stability. It is the 
angular resolution that is most crucial in determining when the numerical instability be- 
comes serious. Higher angular resolution does give more accurate results, but also leads to 
instabilities at earlier times. 

We have performed a series of convergence tests on our code. Convergence was measured 
along the line 6 = tt/4 for a number of metric functions, and for the conformal factor. 
Because we did not have data placed along this value of 6 we interpolated it from our existing 
data using a third order interpolation scheme. The convergence rate of a given quantity was 
calculated by comparing results obtained at three resolutions in a similar manner as reported 
in Ref. [|10|- The basic principle is to assume the error in a given quantity is proportional to 
{Ar]Y, and then the convergence rate a is determined experimentally. For a formally second 
order accurate numerical scheme, such as ours, one expects a ^ 2. 

For the purposes of these tests we required that AO = Ar]. As in previous work |1^ 



we performed most of the tests with a unit lapse. We found, in general, that all quantities 
converged to second order with slight variations throughout the domain. This applied to 
both high amplitude Brill wave data sets and rapidly rotating data sets. In Fig. || we show 
the result of a convergence test for a pure Kerr spacetime with J = 5 (this turns out to have 
a/m = .676). We evolved it at three different resolutions corresponding to 75, 150, and 300 
radial zones and checked the convergence of the radial variable A. The results are shown 
after 4.8M. In two places the direction of the convergence changes, making the convergence 
almost impossible to measure there. These points are each labeled "crossing point" in the 
figure. These results are typical for a variety of convergence tests we have performed. 



V. DISCUSSION 

A. Comparison with ID codes 

In this section we discuss the evolution of the metric functions, both to compare the 
present code to the non-rotating NCSA code, and to show what effect rotation has on the 
system. For these purposes we compare and contrast the evolution of two different black 
hole spacetimes representing a wide range of the kinds of problems our code can evolve. The 
first is a pure Schwarzschild black hole, evolved with a symmetric lapse. This has become a 
classic test problem for black hole codes P, [T0|JT8| . Because the system can be evolved very 



accurately in ID (given sufficient resolution), and also because it is a difficult problem due to 
very large peaks that develop in the solution pSf , it is a strong testbed calculation. In Fig. |^ 
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we show a comparison of our new code for rotating spacetimes with a ID code described in 



Ref. [|TOl- Both codes were run with the same radial resolution, ^ 0.020, with the same 
time steps At = At]. In the 2D case the angular resolution was taken to be AO ~ 0.033. 
The radial metric function A is shown at several times for both codes. Only a single angular 
zone is shown: The initial data are spherically symmetric, and our codes maintain this 
symmetry to a high degree of accuracy throughout the evolution. All other angular zones 
are indistinguishable. The agreement between the two codes is excellent, through the time 
t = lOOM shown here. It is important to note that the 2D evolution was performed as a full 
2D problem without specializing it in any way. For example, the maximal slicing equation 
was solved as a full 2D elliptic equation in the 2D code, whereas it is solved as a simple 
ordinary differential equation in the ID code. 

A key feature of this and the other black hole spacetimes studied here is that constant 
7] observers fall inward toward the singularity, actually passing through the horizon in finite 
time. As a result, more and more grid zones represent regions inside the horizon as the 
evolution progresses. Because of the differential rate at which these observers fall through 
the horizon the grid stretches, creating a sharp peak in the radial metric function A as the 
evolution continues, as shown in Fig. ^. This effect assures us that the code will eventually 
become inaccurate and crash. This is the bane of all black hole evolution evolution codes at 
present, and is the main motivation for considering apparent horizon boundary conditions 
that can, in principle, eliminate this problem |2^j29|] . 



We have also evaluated the error in the Hamiltonian constraint. This error is in the form 
of a mass density which is given by 

P=^{R-KatK''' + K^). (33) 

The error in this quantity is dominated by the axis instability. However, because this insta- 
bility grows most rapidly near the horizon its effect on the code is partially cancelled by the 
collapse of the lapse. Because we are not as interested in violations of the Hamiltonian con- 
straint inside the horizon, since this should not affect what happens outside for hyperbolic 
evolution, we will be more interested in the quantity a\p\. (Although we do solve elliptic 
equations, which propagate information instantaneously, these equations are for gauge con- 
ditions. An "error" in a gauge condition, which is arbitrary, does not affect physical results 
in principle.) The maximum error in this quantity as a function of time for various runs 
is plotted in Fig. 0. Our evolutions are generally unable to proceed past the point where 
this maximum density passes unity. The Hamiltonian constraint violation illustrated by the 
maximum error is dominated by the axis instability, and highly localized. 

Because of the strong locality of this error and the hyperbolic nature of our evolution, 
we do not feel, however, that the maximum error provides us with the best understanding 
of the overall accuracy of the code. For this reason we also consider the average value of 
IpI over the grid (weighted by the lapse function a). We provide a plot of this measure of 
the error for the same runs in Fig. |^. It should be noted that this quantity is generally 4 - 
8 orders of magnitude smaller than the maximum error on the grid for all runs considered. 
The runs that represent distorted spacetimes are initially less smooth, and for this reason 
the Hamiltonian violation, is larger than for Schwarzschild. However, it is small and grows 
at a rate comparable to that of the Schwarzschild evolution. 
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B. Comparison with the 2D codes 



Although the ID problem is a good test of the longitudinal component of the field it does 
not test the code's ability to handle a highly non-spherical, distorted black hole. The next 
case we consider is a distorted, non-rotating black hole with a Brill wave, labeled run rO 
and is specified by (Qo, c", Ji = (1.0, 1.0, 1.0, 0.0, 2). Initially, this data set is a highly 
distorted and nonspherical black hole that evolves in an extremely dynamic way, so it is a 
very strong test case for this new code. In Fig. ^ we show the evolution of this distorted non- 
rotating black hole spacetime obtained with our code, and compare it to the same evolution 
obtained with the code described in Ref. P,pi]|. The radial metric function A = g^^]/'^'^ is 



shown at time t = 60M. Only a single angular zone is plotted (^^ = 7r/2), as by this time the 
developing peak is quite spherical and nearly independent of 6'. In the figure we compare 
the evolution obtained with the rotating code with symmetric lapse, antisymmetric lapse 
conditions, and the non-rotating code of Ref. with symmetric lapse. When the same 

lapse conditions are used, both codes show excellent agreement. Note also that even with an 
antisymmetric lapse (though the early evolution is not shown) the runs will be noticeably 
but minimally different at late times. This is the first time these data sets have been evolved 
with this antisymmetric slicing condition. 

We have compared many other aspects of results of our code with those obtained with 
the code described in Ref. P,p!0[| and find excellent agreement for metric functions, extracted 
waveforms, horizon masses, etc. The comparisons have been performed for a wide range of 
data sets. A detailed study of horizons and waveforms will be presented in Paper II. 



C. Rotating spacetimes 

Rotating spacetimes differ only slightly from the picture discussed above. Constant r] 
observers along the equator are, in general, rotating about the black hole. Because of this, 
their fall through the horizon is slower and grid stretching is less near the equator. The 
peak in the metric function A is, therefore, lesser in magnitude here. This gives the data 
for a rotating black hole spacetime an obvious angular dependence not present in distorted 
non-rotating metric data. 

The pure stationary Kerr spacetime is, in general, much less stable numerically than the 
Schwarzschild, or even the distorted Kerr, and special care is required to evolve it. We have 
two methods which we use to carry out the evolution. The first is to use a lapse which 
is symmetric across the throat. We begin our discussion of rotating black hole evolutions 
by showing a surface plot of metric variable A = g^rj/^^ for a symmetrically-sliced high 
resolution (300 x 30) J = 5 {a/m = .68) Kerr spacetime in Fig. |10|. The error in J as 
computed by Eq.(^) reached a maximum error of 3.6% over the grid at this late time, and 
occurs over the peak in gf^^. This is an important and nontrivial test, as J must be conserved 
during the evolution even though the different metric and extrinsic curvature components 
entering Eq.(^) evolve dramatically. 

Next we consider this same spacetime with an antisymmetric lapse. The antisymmetric, 
maximal sliced Kerr spacetime is actually a time independent analytic solution, but it is 
numerically unstable. Any slight numerical error will destroy the perfect balance between 
the Ricci terms and the derivatives of the lapse in the extrinsic curvature evolution equations. 
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Then the solution will evolve rapidly, behaving much like the symmetrically sliced Kerr 
spacetime. In order to evolve this system stably for long times, we force metric variable F = 
Qecf,/'^^ to be zero at all times (In this stationary spacetime this must be true). Extending 
the grid to larger radii is also helpful. When both of these things are done we are able 
to evolve past 50M before serious problems develop. The J = 5 Kerr evolution with an 
antisymmetric lapse at the throat begins to look like the symmetric lapse run at late times. 
A peak forms in the radial metric function A, although a much smaller one than is seen with 
the symmetric lapse. This is illustrated in Fig. 0. Because the lapse is collapsed in the inner 
region a is virtually zero over a number of grid zones and this satisfies both the symmetric 
and antisymmetric conditions to a good approximation. Note that even the Schwarzschild 
spacetime is difficult to evolve with an antisymmetric lapse since the axis instability becomes 
quite serious at about 60M. The error in J as computed by Eq. (pb]) reached a maximum 
error of 1.6% over the grid. 

Finally, we present a distorted rotating black hole. Calculation labeled run is specified 
by (Qo; ''70; 0"; J) ^) = (1-0, 1-0, 1.0, 10.0, 2). Thcsc labels correspond to the same simulations 
discussed in detail in Paper II, where many physical properties of the spacetimes are ana- 
lyzed. On a 300 x 48 grid, the maximum error in the momentum integral was 1.6%. Again, 
at this high angular resolution the axis instability develops rapidly after about t = TOM, 
causing the code to crash. Lower resolution runs can be carried out past t = lOOM with 
similar results. 

In Fig. |l^ we show a surface plot of the radial metric function A for the rotating hole 
at time t = 60M. Although the familiar peak develops in this function, it does develop the 
expected angular dependence. This is typical of all our rotating black hole evolutions. 

Another typical feature of these spacetimes is "slice-wrapping," discussed in sec- 
tion |1 V A 1| . This name refers to the fact that slices inside the horizon approach a limit 
surface with roughly constant r value. This is connected to the discussion of the limit sur- 
faces above, but here we use this feature to illustrate the full 2D behavior of the metric 
functions to show that they evolve as expected. As a result of slice-wrapping, in our space- 
times the value of "^^y/D (which is the circumferential radius when evaluated at ^ = 7r/2) 
becomes constant over larger regions of the spacetime as the evolution continues. In a non- 
rotating spacetime with a lapse that is symmetric across throat this value goes to 3M/2. 
These spacetimes have a different limiting value, as discussed in section [IV A| . In Fig. ^ we 
show an example of this effect for the spacetime labeled r4. The slice-wrapping effect is clear 
inside the radius of about t] ^ 2.7 where the slice moves out away from the limit surface, 
out across the horizon. From studies of horizons of black holes like those in Ref. pO|, we 
know the horizon is located near t] = 2.7. 

The strong internal structure evident in D = g^^/"^^ well inside the black hole, near the 
throat (?7 = 0), is a remnant of the initial Brill wave, indicating that this black hole was 
initially quite distorted. This structure developed early on as gravitational waves propa- 
gated into the hole, but was "frozen in" as the lapse collapsed rapidly in this region of the 
spacetime. 

A similar behavior is observed in the metric variable B = ggg/'^^ as in except that 
while D = 1 initially, B = e'^^''~'^°^ and so the antisymmetric lapse "freezes in" these different 
functions at the throat. Because the spacetime can evolve even a short distance away from 
there a sharp peak develops at the throat with this slicing condition. Fortunately, the lapse 
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is always small in this region and the code does not suffer as a result. In Fig. we show 
the behavior of B, and Fig. 15 we show the behavior of D. In both plots we are looking 
at run labeled r4 at time 60M. At high resolution (300 x 48) the axis instability becomes 
serious at about TOM, but with lower resolutions it can go past lOOM. 



VI. CONCLUSIONS 

The study of rotating spacetimes presents a new level of complexity to the distorted 
black hole, and we have developed a new code to evolve such spacetimes by building on 
previous non- rotating work [0. Although the rotation requires the introduction of new 
metric and shift variables, we have been able to bring many of the same numerical methods 
to bear on the problem that have been used before. We have shown that the new code 
is able to reproduce results of Ref. [p|,pip!0| for spherical and highly distorted non-rotating 
holes, including the behavior of both the metric functions and derived quantities. We have 
also shown the effect of rotation on the metric functions and how they behave differently for 
rotating black hole spacetimes. In addition to reproducing previous results, the code has also 
passed other tests, such as convergence tests and the conservation of angular momentum, 
and it confirms the theoretical predictions for the relation between the limit slice for a 
maximally sliced black hole and its rotation parameter. 

We have also introduced a new family of distorted black hole data sets, including distorted 
rotating black holes and odd-parity distorted non-rotating black holes. The next paper in 
this series will discuss results obtained from applying this code to these new data sets. 
Specifically it will show how to calculate the location of the apparent horizon and how to 
extract the waveforms for the various i modes radiating from the black hole. It will also 
analyze the behavior of these aspects of the spacetimes for a series of evolutions forming 
a sequence of black holes with increasing rotation, and also for odd-parity distorted black 
holes. 
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FIGURES 



FIG. 1. This plot shows rc, the circumferential radius of the maximally sliced black hole as a 
function of radial coordinate 77, for run labeled ol at time t = TOM. The theoretical limit surface 
for a non-rotating hole is r = 1.5M, but our boundary condition on the lapse prevents us from 
reaching that point at the throat (77 = 0). Nevertheless, away from the throat, the predicted limit 
surface is reached. 

FIG. 2. In this figure the circumferential radius rc is plotted at time t = 60M. For Kerr black 
holes the limit surface is not Tc = 1.5M but a higher value that depends on the rotation parameter 
a/m, as predicted by Duncan. The dashed line is the theoretically predicted r = constant limit. 
In regions away from the throat (where the lapse is zero and the spacetime may not evolve) we see 
that this limit surface is reached. 

FIG. 3. We show a 2D plot of the shift (3'^ as computed by the rotating code for a Kerr black 
hole with J = 5. It is virtually identical to the analytic shift function one would obtain from the 
exact Kerr solution. 

FIG. 4. In both plots that follow, all of the angular zones are displayed, (a) This 2D plot 
shows the contents of metric variable B = gge/^'^ at a resolution of 150 x 24 for run r4., a highly 
distorted rotating black hole. Using this resolution, we can evolve the spacetime to lOOM before 
numerical instabilities develop, (b) This 2D plot shows the metric variable B = gee/^^ (where 
^E* is a conformal factor) in run r4- It has a higher spatial resolution (300 x 48) and it develops 
trouble at an earlier time. The ridges near the throat are a result of the axis instability and the 
radial integration we use to compute the (i'^ component of the shift. 

FIG. 5. This figure plots the convergence exponent a as a, function of radial coordinate r] 
for the radial metric variable A = g^y^/^^ (where ^' is a conformal factor) of a Kerr black hole 
spacetime with J = 5, as measured along 9 = tt/4, evolved to t = 4.8M. As discussed in the text, 

a was measured by evolving the spacetime at 3 different solutions. As our methods are second 
order, we expect to see roughly a horizontal line at cr = 2 (cr is the order of convergence). The 
points labeled "crossing point" are places at which the direction of the convergence is changing, 
and we do not expect to be able to measure the convergence well at these points. 

FIG. 6. This plot compares a ID code to our 2D code for a Schwarzschild spacetime. Each solid 

(dashed) line represents the metric variable A = g^^^/'^'^ (where r/ is a logarithmic radial coordinate 
and ^' is a conformal factor) from the 2D (ID) code at a time interval of lOM, beginning with 
lOM and continuing to lOOM. 

FIG. 7. In this plot we show the maximum violation in the Hamiltonian constraint function 
a.\p\ on the grid. The error is dominated by the effects of the axis instability and is highly local- 
ized there, and evolution becomes impossible shortly after this error exceeds unity. The error is 
measured in units of M^^j^. 
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FIG. 8. In this plot we show the average of the Hamiltonian constraint function \p\ (weighted 
by the lapse function) . The errors for distorted black holes is noticeably larger since these functions 
are initially less smooth. The error is measured in units of M^^j^. 

FIG. 9. In this figure we compare several calculations of the radial metric function A at time 
t = 60M for run rO (a distorted black hole with no angular momentum nor odd-parity distortion). 
Three results are shown: both codes were run with a symmetric lapse across the throat, and the 
present rotating code was also run with an antisymmetric lapse across the throat. We see that 
the rotating and non-rotating code produce practically identical results for the symmetric lapse. 
We also see that the antisymmetric lapse produces similar results, except near the throat where it 
must be different. 

FIG. 10. This surface plot shows the radial metric function A = g^p^/'^^ (where r) is a loga- 
rithmic radial coordinate and ^' is a conformal factor) at time t = 80M for an antisymmetrically 
sliced Kerr spacetimc with J = 5 (a/m = .68). The characteristic peak in A is well-developed and 
its angular dependence, explained in the text, is clearly visible. 

FIG. 11. This surface plot shows the radial metric function A = g^^i/^'^ (where 77 is a logarith- 
mic radial coordinate and is a conformal factor) at time t = 50M for an antisymmetrically sliced 
Kerr spacetime with J = 5 (a/m = .68). At this high resolution (300 x 30), the axis instability 
sets in shortly after this time. 

FIG. 12. This is a surface plot of metric function A = 5^^/^^ (where is a logarithmic radial 
coordinate and ^' is a conformal factor) at time t = 60M for the distorted rotating black hole 
run labeled r4- The presence of a shift in the evolution means rj = constant observers near the 
equator feel a "centripetal force" and fall in more slowly. This results in less grid-stretching, and 
a lower peak there. 

FIG. 13. This surface plot is of the function R = ^r^D^/^ for the run labeled r4. 
Along the equator this value is the equatorial circumferential radius. This plot shows that the 
"slice-wrapping" effect affects a large portion of the grid. 

FIG. 14. This surface plot depicts metric variable B = gee/'^* (where ^' is a conformal factor) 
at time TOM for the run labeled r^. Note that the angular shape of the Brill wave on the initial 
data slice is preserved at the throat, but outside this region it becomes very close in value to metric 
variable D = g^^/^>'^ . 

FIG. 15. This surface plot depicts metric variable D = 5<^<^/\I'"^ (where ^' is a conformal factor) 
at time 60M for the run labeled r4- Note that at the throat D = g^pfj)/'^^ has value 1, as it does in 
the initial data. Farther from the throat it has the same value as B = goo/'^'^. 
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